Skewed Multivariate Bimbaum-Saunders Distributions 

Artur J. Lemonte 

Departamento de Estatistica, Universidade de Sao Paulo, Sao Paulo/SP, Brazil 

Guillermo Martmez-Florez 

Departamento de Matemdticas, Universidad de Cordoba, Monten'a, Colombia 

German Moreno- Arenas 

Escuela de Matemdticas, Universidad Industrial de Santander, Bucaramanga, Colombia 

Abstract 

The univariate Birnbaum-Saunders distribution has been used quite effectively to model times 
to failure for materials subject to fatigue and for modeling lifetime data. In this article, we define 
a skewed version of the Bimbaum-Saunders distribution in the multivariate setting and derive 
several of its properties. The proposed skewed multivariate model is an absolutely continuous 
distribution whose marginals are univariate Bimbaum-Saunders distributions. Estimation of the 
parameters by maximum likelihood is discussed and the Fisher's information matrix is deter- 
mined. A skewed bivariate version for the generalized Bimbaum-Saunders distribution is also 
introduced. We provide an application to real data which illustrates the usefulness of the proposed 
multivariate model. 

Key words: Birnbaum-Saunders distribution, generalized Bimbaum-Saunders distribution, max- 
imum likelihood estimators, modified moment estimators, multivariate distributions. 



1 Introduction 



The univariate family of distributions proposed by iBimbaum and SaundersI (|1969|) . also known as 
the fatigue life distributions, has been widely applied for describing fatigue lifetimes. This family 
was originally derived from a model for which failure follows from the development and growth of a 
dominant crack. A random variable T has a Bimbaum-Saunders (BS) distribution if it can be written 
as T = f3{aZ/2 + + 1]^/^}^, where Z is a random variable following the standard normal 

distribution, i.e. Z ~ N(0, 1). Its density function is 
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which depends on two parameters: the shape a > and scale /3 > 0, which is also the median of the 
distribution. We have kT ~ BS(a, k/3) for any A; > 0, i.e. the BS distribution is closed under scale 
transformations. The expected value, variance, skewness and kurtosis of T are, respectively. 
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The density function ([T]) is right skewed and the skewness decreases with a. Notice that both mean 
and variance increase as a increases. It is also of interest to mention that if T ~ BS(a,/3), then 
~ BS{a, /3~^). It impl ies that the BS d istribution also belongs to the family of random variables 



closed under reciprocation ([Saunders 



19741) . It then follows that 
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The shape of the hazard function of the BS distribution is discussed in 
authors showed that the hazard rate function is not monotone and is unimodal for all ranges of the 
parameter values. Some interest i ng res ults on improved statistical infer ence for the BS distribution 
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may be revised in lWu and WongI (|2004|) and 

The univariate BS distribution has received signifi cant attention over the l a st few years by maiiy 
researchers and some generalizations are proposed in Diaz -Garcia and Leival (2005 ). Owen ( 2006 ). 



Guiraud et al. 



mm. 
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Castillo et all (2011b and 



Cordeiro and Lemonte (2011 



among other. On the other hand, as far as we know, little wor k has been done to extend the BS dis 
tributi o n to the multivariate case . We can refer to th e work s by 



(2006), 



Kundu et al 



(1201 oh and lCaro-Lopera et all (120 12h . In 



Dfaz-Garcfa and Dommeuez-Molina 



Dfaz-Garcfa and Domfnguez-Molina 



(I2006h . the authors defin ed an independent m ultivariate BS distribution. By using the bivariate nor- 
mal distribution function, Kundu et al. ( 20101) proposed a bivaria t e BS d istribution which is absolutely 
continuous and has five parameters. Finally. ICaro-Lopera et al.l (|2012|) introduced the matrix- variate 
generalized BS distribution. 

As can be observed, little work on multivariate versions for the BS distribution have been pub- 
lished. In this paper, in addition to the existing multivariate BS mo dels, we shall propos e the asym- 
metric (skewed) multivariate BS distribution based on the work of Arnold et al.l (2002). The main 
motivation for introducing this multivariate version of the BS distribution relies on the fact that the 
practitioners will have a new multivariate BS model to use in multivariate settings, since the formu- 
lae related with the new multivariate model are manageable and with the use of modem computer 
resources and its numerical capabilities, the proposed model may prove to be an useful addition to the 
arsenal of applied statisticians. Additionally, the new model is quite flexible (see Figure [Din Section 
|2]) and can be widely applied in analyzing multivariate data. Further, we provide an application to real 
data in which is showed that the new multivariate model yields a better fit than other multivariate BS 
distributions available in the literature. 
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The paper unfolds as follows. The skewed bivariate BS distribution is defined in Section[2]and then 
several properties are discussed. The multivariate extension is presented in Section|3] In Section|4l we 
propose different methods for estimating the unknown parameters as well as derive the information 
matrix and discuss likelihood ratio tests for some hypotheses of interest. In particular, we propose 
modified moment estimators for the unknown parameters which are explicit in form and can therefore 
be used effectively as the initial guess in the iterative process for the computation of the maximum 
likelihood estimators. Further, the asymptotic distribution of the maximum likelihood estimators is 
derived and thus the asymptotic confidence intervals for the unknown parameters can be constructed. 
The usefulness of the proposed model is illustrated in an application to real data in Section|5l We also 
introduce in Section |6]the skewed bivariate generalized BS distribution. Finally, Section |7] closes the 
paper with some concluding remarks. 



2 Skewed bivariate BS distribution 

We initially consider the skewed bivariate BS distribution. For each a; G M and for each y G M, 
consider the conditional distributions 

X|r = y ~ SN(A?/), y|X = X ~ SN(Ax), (2) 

wher e X\Y = V y SN(A?/) means that given Y = y, X\Y = y has skew normal distribu- 



Amold et al. 



tion (I Azzalini ,119 8 51) . The shape parameter A G M determines the skewness of the density. From 
(|2002l) and using the conditional distributions in the joint probability density func- 



tion (pdf) of the random vector {X, Y) takes the form 

/x,y(x, y) = 2(Pix)(l>{y)^Xxy), (x, y) G M^ (3) 

where and $(■) denote the pdf and cumulative distribution function (cdf) of the standard normal 
distribution, respectively. Also, fx{x) = (l){x) and fyiv) = 4>{y)- If A = in Q, then fxxi^^v) = 
(j){x)(j){y) and hence X and Y become independent. For A 7^ 0, it can be shown that the correlation 
between X and Y, p{X, Y) say, is given by 

(^v^ ■ (x^ ^(3/2,2,l/(2A^)) 
piX,Y) = sign(A) X , 

where U{a, b, z) denotes the confluent hypergeometric function, defined as 

1 /"°° 

t/(a, b,z) = -— e-^Y-\l + tf-'^-'dt, 

r(a) Jo 

with 6 > a > and z > 0, and r(-) represents the gamma function. Therefore, the parameter A also 
governs the correlation. 
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Let Zj ~ N(0, 1), for J = 1,2, with ZijZa = ^2 ~ SN(A^2) and Z2IZ1 = zi ~ SN(A2i). Then, 
taking the transformation 



T, = /3, 



J = 1,2, 



where aj > and (3j > 0, the joint pdf of the skewed bivariate BS (SBVBS) distribution takes the 
form 



fn,T2{ti,t2) = 20(ai)0(a2)$(Aaia2 



where 
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The notation used is (Ti, T2) ~ SBVBS(ai, ^2, /32, A). The random variables Ti and T2 become 
independent for A = in dH) and hence the proposed bivariate model redu ces to the independent 
bivariate model considered by Diaz-Garcia and Dommguez-Molinal (|2006h . So, as remarked, the 
shape parameter A also introduces correlation between Ti and T2. 

Contour plots for the joint pdf dU are presented in Figure [T] From this figure, note that dU can 
take on different shapes and will therefore be useful in analyzing bivariate data. Additionally, notice 
that @) can be unimodal or bimodal depending on the value of A. 

The following theorem provides the marginal and conditional distributions of the SBVBS distri- 
bution. 

Theorem 2.1. //(Ti, T2) ~ SBVBS(ai, ^2, /32, A), then: 

(i) T, ~BS(a,-,/3,),/orj = l,2. 

(ii) The conditional pdfofTi given T2 = ^2 is 



fn\T2{ti\T2 = t2) = 20(ai)$(Aaia2; 



(iii) The cdfofTi given T2 = ^2 is 



Pr(Ti < ti\T2 = t2) = $(ai) - 2T(ai, Aa2), 



where T(-, ■) denotes the Owen's function aOwen . 



1956 ). 



Proof. Parts (i) and (ii) follow from de definition of the distribution. We have that 

Pr(Ti <ti|T2 = t2) = 



ti J.-3/2/J- \ fj ) 
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(b) 








(f) 
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Figure 1: Contour plots of the density function dU for some values of (ai, 02, A, /32, A): 
(a) (0.5,0.5,1.0,1.0,0.5); (b) (0.2,0.4,1.0,1.0,-1); (c) (0.8,0.5,1.0,1.0,1.5); (d) 
(1.5, 1.5, 1.0, 1.0, 1.5); (e) (0.2, 0.2, 1.0, 1.0, 5); (f) (0.2, 0.4, 1.0, 1.0, -10); 
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where 



at = — 



1/2 



t 



1/2' 



Making the change of variable u = at,we arrive at 

Pr(Ti < ti|T2 = ^2) = / 20(M)$(AMa2)dM. 



Now, from 



Azzalinil (119851) we can show that Pr(Ti < ti|T2 = h) = $(ai) - 2T(ai,Aa2) and 
therefore the result (iii) holds. □ 

Some properties of the random vector (Ti, T2) are provided in the following theorem. 

Theorem 2.2. //(Ti, T2) ~ SBVBS(ai, ^2, /3i, /32, A), 

(i) (fciTi, r2) ~ SBVBS(ai, a2, /32, A), fci > 0. 

(ii) (Ti, A;2T2) ~ SBVBS(«i, ^2, A, A;2/32, A), A;2 > 0. 

(iii) (fciTi, A;2r2) ~ SBVBS(«i, ^2, A;2/32, A), h, h > 0. 

(iv) {T-\ T,') ~ SBVBS(«i, «2, /3r\ /32'\ A). 
(V) (Tr\T2) ~SBVBS(«i,«2,/3r\/32,-A). 
(vi) (Ti, T^-i) ~ SBVBS(«i, «2, /3i, /32"\ -A). 

Proof. These follow from dH) upon using suitable transformations. □ 

Since the marginal distributions of the bivariate vector (Ti, T2) are BS distributions, the mean and 
variance of Ti and T2 are obtained directly from these marginals in the forms 



E(T,) = /3, (1 + , V(T,) = (l + , 



J = 1,2. 



Additionally, 



J = 1,2. 



The product moments of (Ti, T2), E(T{'T|) say, are very complicated to be determined algebraically 
and have to be computed numerically. In the following, we shall derive an expression for E(TiT2) 
which can be of some interest. We can show after some algebra that 
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being ci = 1, C2 = 4A^/ (3!) and cs = 32A^/ (5!). For A = (independent case), we have immediately 
that 



3 Multivariate extension 



We have considered the bivariate case in Section |2l but extensions to higher dimension can be readily 
accomplished using suitable notation. For a random variable Z = (Zi, . . . , ZpY of dimension p, we 
define the subvectors -^(i),. • • ,Z{p) of dimensions (p— 1) such that, for each j = 1, . . . ,p, Z(^j) denotes 
the vector Z with the jth coordinate Zj deleted. Analogously, for a real vector z = (zi, . . . , Zp)~^ , 
Z(^j) is obtained from z with the jth coordinate zj deleted. 
By assuming (for each j = 1, . . . ,p) that 



the joint pdf of Z = (Zi, . . . , Zp)^ takes the form (jAmold et al 



2002) 









fz{z) = 2 
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Thus, under the transformation 

T, = /3, 
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where Zj ~ N(0, 1), we obtain the joint pdf of T = (Ti, . . . , T^Y in the form 



n ^^^i' 



^n^. n 

j=i / j=i 



where > 0, /3j > and is given in Q, j = 1, . . . ,p. Let ck 



t e 



(6) 



and /3 



. . . , If T = (Ti, . . . , Tp)^ has skewed multivariate BS distribution, then we use the notation 
T~ SMVBS(a,AA). 

Several properties discussed in the bivariate case hold for this multivariate extension. For exam- 
ple, Tj ~ BS(aj, f3j) for j = 1, . . . ,p, i.e. the marginal distributions are BS distributions; A = 
corresponds to the independent case; for ki, kp > 0, (fciTi, . . . , kpTp) ~ SMVBS(q:, /3*, A) with 
f3* = . . . , kpPpy- (Tf 1, . . . , T-i) ~ SMVBS(a, f3**, A), where f3** = (/JfS • • • , f3;'V, and 

so on. In the next section, we shall consider estimation for the unknown parameters of the SMVBS 
distribution in as well as inference. Thus, from these general results the bivariate case considered 
in Section [2] can be easily specialized by considering p = 2. 



4 Estimation and inference 

In this section, we address the problem of estimating the unknown parameters of the SMVBS dis- 
tribution. Let ti, . . . ,tn denote a random sample of the SMVBS(q;, /3, A) distribution, where ti = 
(tii, . . . , tpiy and, as before, a = (ai, . . . , ap)^ and (3 = . . . , /3p)^. Let 9 = {a.^ , f3^ , A)""" be 
the parameter vector of interest of dimension 2p + 1. 



4.1 Modified moment estimators 



First, we shall prese nt modified mom ent estimators (MMEs) for the unknown parameters by follow- 
ing the approach of Ng et al.l (l2003h . The SMVBS model has 2p + 1 parameters and the marginal 
distributions are BS distributions with parameters (aj, (3j), j = 1, . . . ,p. Then, the moment estima- 
tors for aj and can be obtained by equating ]E(Tj) and V(Tj) to the corresponding sample estimates 
for J = 1, . . . ,p. However, it is known that in the case of univariate BS distribution, the moment esti- 
mators may not always exist ( Ng et al. . 12003 ). Here, we will use lE(Tj) and E(Tj~^) instead of using 
E(Tj) and V(Tj), and equate them to the corresponding sample quantities. After some algebra, the 
MMEs for ai, . . . , Op and . . . , /3p are 



1/2 



1/2 



4- 



,p, 



(7) 
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The MMEs for cx and /3 in (|7]) are explicit in form and can be used effectively as the initial guess in 
the iterative process for the computation of the maximum likelihood estimators (MLEs) in the next 
section. 



4.2 Maximum likelihood estimators 



The log-likelihood function for the parameter vector (apart from an unimportant constant) is given 
by 



i{e) = -nJ2 

^ n p 



n p 



i=l j=l 



1=1 



i=l j=l 

p 



(8) 



where 



1/2 1/2- 
tji 



The MLEs of the unknown parameters are obtained by maximizing the log-likelihood function in ^ 
with respect to 6. By taking the partial derivatives of the log-likelihood function in ([8]) with respect 
to the parameters aj, /Sj and A, we have (for j = 1, . . . , p) 



di{e) 



n 1 " A " ^ 
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where 



Wi = Wi{a.,(3,X) 



di{e) 
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i=i j=i 
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TheMLE^ = A)^ of6 = {a^ ,/3^ , A)^ can be obtained by solving the likelihood equations 

d£{e) d£{0) de{e) 



0, j = l,...,p, 



(9) 
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simultaneously. There are no closed form expressions for the MLE and its computation has to be per- 
formed numerically using a nonlinear optimization algorithm. The Newton-Raphson iterative tech- 
nique could be applied to solve the likelihood equations and obtain the estimate 6. For computing 
the MLEs, starting values for the algorithm are required. Since the MMEs for aj and /3j in (|7]) are 



explicit, they can be use d effectiyely as the initial guess in t he iterative procedure. The 0^ rn atrix 



programming language (IDoomikl 120061) and the R program (|R Development Core Team , 
be used to compute 6 numerically. 

We can show from the likelihood equations that, for given I3i, . . . , fip, the MLEs of ai, 



2010) can 



Op are 



- 2 



1/2 



,P- 



By replacing aj by aj{j3j) in ([8]), we obtain the profile log-likelihood function for /3 and A as 



log(a,(/3,)) + -log(/3,) 



ep{f3,X) = -nJ2 

^ n p 



n p 



i=i j=i 



i=i j=i 



2=1 



where 



1/2 



1/2 



We can also obtain the MLEs of (3 and A by maximizing the profile log-likelihood function A) 
with respect to /3 and A. The Newton-Raphson algorithm or some other optimization algorithm to 
maximize £p(/3, A) with respect to /3 and A needs to be used, since the MLEs of /3 and A cannot be 
obtained explicitly. The profile log-likelihood function ^p(/3, A) is not a real log-likelihood function 
and some of the properties that hold for a genuine log-likelihood do not hold for its profiled version. 
In particular, there exist score and information biases, both of order 0(1). 

The asymptotic inference for the parameter vector 6 = (a^, A)^ can be based on the nor- 
mal approximation of th e MLE 6 of 6 = (ck^, /3^, A)^. Under some regular conditions stated in 
Cox and Hinkleyl (| 19741 . Ch. 9) that are fulfilled for the parameters in the interior of the parameter 
space, we have 6 ~ N2p+i(^, S^^), for n large, where ~ means approximately distributed and 
is the asymptotic variance-covariance matrix of 0. The matrix Eg is given in the Appendix. The 
multivariate normal N2p+i(0, Sg^) distribution can be used to construct approximate confidence in- 
tervals for the parameters aj, (5j and A, which are given, respectively, by cij ± x [¥(5?^)]^/^, 
± z^i2 X [V(;5j)]i/2 and A ± z^i2 X ^(X)^/'^, where V(-) is the diagonal element of ^ available 
at 6 corresponding to each parameter, and z^j2 is the quantile 100(1 — 7/2)% of the standard normal 
distribution. 



'Ox is freely distributed for academic purposes at ^ttp://www.doornik.com 
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Besides estimation of the model parameters, hypotheses tests can be taken into account. Let 
d = (6j , Oj)^ , where 0i and 02 are disjoint subsets of 6. Consider the test of the null hypothesis 
Ho '■ 01 = 0m against 1-Li : 0i 0qi, where is a specified vector. Let be the restricted MLE of 
obtained under "Ho- The likelihood ratio (LR) statistic to test "Ho is given hy u = 2{£(0) — (■{0)}- 
Under "Ho and some regularity conditions, the LR statistic converges in distribution to a chi-square 
distribution with dim(0i) degrees of freedom. In particular, the LR statistic to test the null hypothesis 
7^0 : A = against "Hi : A 7^ takes the form 

u; = 2{£(a,AA)-£(5,A0)}, 

where a. and (3 are the restricted MLEs of ol and (3, respectively, obtained from the maximization 
of dHl) under I-Lq : \ = The limiting distribution of this statistic is Xi under the null hypothesis. 
The null hypothesis is rejected if the test statistic exceeds the upper 100(1 — 7)% quantile of the Xi 
distribution. 



5 Application to real data 



In this section, for illustrative purposes, we present an empirical application to demonstrate the ap 
plicability of the proposed skewed m ultivariate BS distribution. For the sake of com par ison, we also 



consider the distributions proposed in iDiaz-Garcia a nd D o mmguez-Molinal (|2006|) and 



(|2010l) . We shall use the data set obtained from 



Kundu et al 



VoUa (|1985|) . which represent the amount of time (in 



hours) spent on two categories of activities over 100 days in the year 1976 for 28 individuals. The 
data are: (115, 175), (100, 115), (130, 160), (115, 180), (119, 143), (100, 150), (960, 132), (150, 115), 
(142, 870), (180, 125), (152, 122), (174, 119), (140, 100), (147, 840), (105, 700), (950, 600), (130, 
600), (105, 800), (117, 650), (850, 400), (102, 450), (100, 960), (920, 640), (128, 860), (102, 122), 
(107, 730), (860, 580), (940, 580). The first figure represents the amount of time spent on eating and 
the second figure represents the amount of time spent on watching tele vision. All the computations 



were done using the Ox matrix programming language (iDoomikL I2006|) . 

We now use the SBVBS distribution to model these bivariate data. We obtain from the data 
si = 118.14, S2 = 99.43, n = 113.40 and fa = 84.61, and hence the MMEs are cti = 0.2035, 
02 = 0.4099, $1 = 115.7457 and /32 = 91.7220. These values are used as initial guesses for ai, 
02, (3i and (32, respectively. An initial guess for A is also required to start the maximization of the 
log-likelihood function ([8]), i.e. to solve the likelihood equations dH) with p = 2. As initial value for 
A we consider A = 0, which corresponds to the independent case. The algorithm converges after 
21 steps and the MLEs of ai, 02, A, /32 and A are Si = 0.2047, 02 = 0.4101, /3i = 113.2907, 
^2 = 90.7447 and A = 0.8806, respectively. Notice that the MMEs for ai, 02, /3i and /32 are close 
to their respective MLEs. We have also considered other initial guesses for A, for example, with the 
initial values A = —5,-2,3 and 4, the algorithm converges to the same estimates after 39, 26, 25 
and 31 steps, respectively. The 95% asymptotic confidence intervals for ai, a2, (32 and A are 
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(0.1508, 0.2586), (0.3051, 0.5152), (104.8325, 121.7489), (77.8975, 103.5919) and (0.0349, 1.7263), 
respectively. 

Next, we make use of the LR statistic to test the null hypothesis T-Lq : X = against "Hi : A 7^ 0. 
Here, uj = 2{i{ai, a2, /32, A) — i{ai, 02, /32, 0)}, where 5i, 52, /3i and /32 are, respectively, the 
restricted MLEs of ai, a2, /3i and /32 obtained under "Hq and are given by Si = 0.2035, 0.2 = 0.4099, 
Pi = 115.7470 and ^2 = 91.7128. By a little computation, we have that the LR test statistic (uj) 
equals 6.6834 (p-value < 0.01). Therefore, the null hypothesis T-Lq : A = is strongly rejected at 
the usual significance levels and hence the assumption of the skewness (correlation) is suitable for 
the cu rrent bivariate data. Since the bivariate distribution in iDfaz-Garcfa and Dominguez-Molina 



(|2006|) . DG-DM say, and our proposed model are nested models (i.e. the DG-DM model holds for 
A = 0), the null and altemative hypotheses can be rewritten as T-Lq: DG-DM against Hi: SBVBS. 
Thus, based on the LR statistic above, the SBVBS distribution fits the data better than the bivariate 
DG-DM model. 

The generalized LR statistic {Tlr^nj^) presented in IVuongI (Il989h can be used for discriminat- 
ing among non-nested models, which is a distance between the two models measured in terms of 
the KuUback-Lieb ler information criterion. Then, our proposed model and the bivariate model in 
(120 101) can be compared by using Tlji^^n. For strictly nonnested models, T^r^nn con- 



Kundu et al 



verges in distribution to a standard normal distribution under the null hypothesis of equivalence of 
the models and the null hypothesis is not rejected if |Tiij,ArAr| < z^/2, where z.y/2 is the quantile 
100(1 — 7/2)% of the standard normal distribution. On the other hand, we reject at significance level 
7 the null hypothesis of equivalence of the models in favor of the SBVBS model being better (or 

Kundu et al. i2010 ) if T^r^nn > (or Tlr,nn < —z-X The generalized 



worse) than the model in 

LR test statistic (Tlr,nn) equ als 4.0903 (p-yalue < 0.01). Therefore, the proposed model is signifi- 

(|2010h according to the generalized LR statistic to model 



Kundu et al 



cantly better than the model in 
the current data. 

A natural question at this point is whether SBVBS model fits the current data satisfactorily. Here, 
in order to verify it, we computed the modified Cramer- von Mises (W*) and Anderson-Darling (A*) 
statistics for the fitted marginals, i.e. BS('0.2047, 113.2907) and BSfO .4101, 90.7447). The statistics 



W* and A* are described in details by 



Chen and Balakrishnan 



1995 ). The values of these statistics 



are 0.0971 (p-value > 0.1) and 0.5680 (p-value > 0.1), and 0.0513 (p-value > 0.1) and 0.3145 (p- 
value > 0.1), respectively. Therefore, based on the marginals, we have that the SBVBS distribution 
can be used effectively in this case. Although it does not guarantee that the bivariate real data will 
have SBVBS distribution, at least it gives an indication that the SBVBS model may be used to analyze 
this bivariate data. 
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6 Skewed bivariate generalized BS distribution 



The univariate generalized BS (GBS) distribution was proposed in iDfaz-Garcfa and Leival (|2005|) . 
which is a highly flexible lifetime model that admits different degrees of kurtosis and asymmetry 
and possesses unimodality and bimodality. The GBS distribution is related to standard symmetrical 
di stributions in R, al so kn own as elliptically conto ured univariate distributions. The reader is referred 



to 



Fang et al. 



(ll99ol) and 



Gupta and Vargal (Il993h for more details about symmetrical distributions. 



For the univariate case, elliptical distributions correspond to all the symmetric distributions in M. 
Specifically, a random variable X has an elliptical distribution if its probability density function is 
given by fx{x) = cg{[x — x G M, where /i G M is a location parameter and > is a scale 

parameter. The function (7 : M — )• [0, oo) corresponds to the kernel of the density of X and c is the 
normalization constant such that fx{x) is a density. The function g{-) is typically known as density 
generator. We then write X ~ E(/i, 0^; gf). 

The notation Z ~ E(0, 1; (7) or Z ~ is used for a random variable Z that follows a standard 
elliptical distribution in M. The pdf and cdf of Z are denoted by /(•) and F{-), respectively, where 
f{z) = cg{z^) and F{z) = f{z)dz. The density generator of the normal, Cauchy, Student-t, 
generalized Student-t, type I logistic, type II logistic and power exponential are, respectively, given 



(27r)-i/2 exp{-u/2), g{u) = {ti{1 + u)Y\ g{u) 



-1/ 



u 



by g{u) 

where > and is the beta function, g{u) = s''/^i?(l/2, r/2)"^(s + (s,r > 

0), g{u) = ce~"(l + e~")"^, where c ~ 1.484300029 is the normalizing constant obtained from 
J^u-^/^g{u)du = 1, g{u) = e-^{l + e'^)-'^ and g{u) = c(A:) exp(-|Mi/(i+fc)), -1 < k < I, 
where c{k) = 1(1 + {k + l)/2)2^+^'^+^'^/\ 

In the following, we shall introduce the skewed bivariate GBS (SBVGBS) distribution. A random 
variable Y follows a standard skew-elliptical distribution in M if its pdf takes the form 



friy) = 2f{y)F{\y), 



y e 



(10) 



We use the notation Y ~ SE(A; (?). If A = in (iTOl) . then the standard elliptical distribution holds, 
i.e. Y ~ E{g). Now, let Zj ~ E{g), for j = 1, 2, with Zi\Z2 = Z2 ^ SE{\z2; g) and ^2!^! = -^i ~ 
SE(Azi; g). Additionally, consider the transformation 

-\ 2 



^Z, 



2 ^ 



+ 1 



J = 1,2, 



where aj > and (3j > 0. Then, from the above transformation and using results due to 
h002) . the joint pdf of the SBVGBS distribution is given by 



Amold et al. 



3/2/ 



{tl,t2 



(11) 



2al^/f3l 2a2V/32 

where aj is defined in ©. If (Ti, T2) follows the SBVGBS distribution, the notation used is (Ti, T2) ~ 
SBVGBS(ai, ^2, 132, A; g). Notice that the joint pdf dH) is a special case of (fTTI) . All extra param- 
eters are considered as known or fixed in (fTH . For example, the degrees of freedom for the Student-i 
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model. The main motivation for this generalization of the SBVBS model presented in Section [2] is 
based on the search for bivariate distributions that are more flexible than the SBVBS model in ana- 
lyzing bivariate data. 

Some properties for this bivariate class of distributions are presented in the following theorem. 
Theorem 6.1. //(Ti, T2) ~ SBVGBS(ai, 03, A, (^2, A; g), then: 

(i) T, ~GBS(«,-,/3,;(7),/orj = l,2. 

(ii) (A;iTi, T2) ~ SBVGBS(ai, ^2, hf3i, /^a, A; g), h > 0. 

(iii) (Ti, ^2^2) ~ SBVGBS(ai, ^2, A, ^2, A; g), h > 0. 

(iv) (fciTi, A;2T2) ~ SBVGBS(ai, ^2, hP^ ^(^2. A; g), h,k2> 0. 
(V) (Tr\T2-i) ~SBVGBS(ai,a2,/5r\/32"\A;(?). 

(vi) (rr\T2) ~SBVGBS(ai,a2,/3r\/32,-A;(?). 

(vii) (Ti,T2-i) ~SBVGBS(ai,a2,/3i,/32"\-A;^7). 

Proof. Using suitable transformations in (fTTI ). these results follow. □ 

From (fTTT ). several news SBVGBS distributions can be obtained. For example, the joint pdf of the 
skewed bivariate BS Student-t model takes the form 



. 1 



2' 2 



where z/j is the degrees of freedom, = (Aaia2)^/[(Aaia2)^ + z^j] and I^(r,s) is the incomplete 
beta ration function. The skewed bivariate BS Cauchy distribution is a special case of the joint pdf 
above when z/i = z/2 = 1. It is evident that other bivariate models can be obtained as, for example, the 
skewed bivariate BS type I (type II) logistic model, skewed bivariate BS power exponential model, and 
so on. Further, extensions to higher dimension can be derived and MLE of the unknown parameters 
can also be considered. These problems can be developed in a future research. 



7 Concluding remarks 

The univariate BS model has many attractive properties and has found several applic ations in the 



literature in cluding lifetime, survival and environmental data analysis (see, for example, iLeiva et al.L 



2008 



20091) . As mentioned before, little work has been done to extend the BS model for the mul- 
tivariate case. In this article, we have introduced the skewed multivariate BS distribution. The new 
distribution is very general, quite flexible and widely applicable. The new model is an absolutely con- 
tinuous multivariate distribution whose marginals are univariate BS distributions. We have discussed 
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several properties of this new class of distributions and the estimation of parameters is approached 
by the method of maximum likelihood. The observed and expected information matrices are deter- 
mined and likelihood ratio tests for some hypotheses of interest are also considered. The skewed 
bivariate BS distribution is discussed and we have shown that the additional shape parameter (A) in- 
troduces skewness, correlation and bimodality to this distribution. These interesting properties make 
this bivariate model a quite flexible d istribution to model bivariate data. Other biv ariat e BS models 



Kunduetal 



have b een introduced and are given in iDiaz-Garcia and Dominguez-Molinal (|2006l) and 
hold) . KB J say. The DG-DM model is an independent bivariate model and hence does not consider 
correlation between the random bivariate vector. The KBJ model considers correlation between the 
random bivariate vector, but does not allow bimodality. As remarked, the skewed bivariate BS model 
proposed in this article can be skewed, correlated and bimodal, and therefore is much more flexible 
than the other bivariate BS models available in the literature for analyzing bivariate data. This is sup- 
ported in an application to real data in which we show that the skewed bivariate BS model provides 
consistently better fit than the DG-DM and KBJ models. Finally, we have also introduced in this pa- 
per the skewed bivariate generalized BS distribution and discussed some of its properties. Although 
we have discussed the generalized BS distribution in bivariate settings, the skewed multivariate gen- 
eralized BS distribution can be introduced along the same lines. This problem can be developed in a 
future research. 
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Appendix. Fisher information matrix 



We present the elements of the Fisher information matrix S^. First, we shall compute the elements 
of the Hessian matrix 



with 



docdoL 



dH{e) 







tT 


Lj3j3 Ljsx 


J-T 


L^x Lxx _ 


L 





■laX 



dH{0) 
dctdX 



(LaiX, ■ ■ ■ , LapX) 



dHio) 



-11313 



d/3df3^ 



■'I3X 



dH{e) 

df3dX 



{Lp^x, ■ ■ ■ ,L(s^xV, 
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where = l,...,p. 



n 3 " 2A " ^ 

j j i=l 3 i=\ j=l 



n p 



- ^ E n 4 - E n ^. 



J i=l j=l '^j 1=1 j=l 



2 



n p 



X 



2 n p 



La^,aj = E n E n 

■' ■' 1=1 j=l i=l 

^3 n p 



1=1 1=1 



1 " 

~ Q,3«. E 
"jPj i=l 

A" 



A 



i=i 



n 



Ojj'ii 



A 



-1 + A^ ]^ 4 + -^^^ n '^i* 



n "j'^' ^ 



^jPj i=i ^Pj i=i j=i 



m .=1 



A 



A 



-1 + - — dij \Wi + XY[ciji]Yl ^fi 

^"^^ ^ 3=1 ) 3'+3 . 



n 



dj'i, 



A 

-- — T 3- Widijdif 



■1 + x[wi + xyIojAyi 



n 



f + J, 



Wi 



1 " 

--E 

i=l 



1 " 



1 - A I A JJ^ + Wj JJ^ 

j=i 3=1 



3=i 



1-Xlxl[a% + Wil[a 



^•^A = - E n 4 + ^» n ) n 

i=i \ j=i j=i J j=i 



n 



Clj'ii 
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The Fisher information matrix is given by 



Hff — —'E,{Lffff) 



^aot ^aX 
^/9/9 ^/3A 



where 



■lQ.13 



for j,/ = 



_ 2n A-^ ^ 

J J 1=1 



n 



IT 



A 



\2 



^'114 



j" J, 



A' 



p p 



O-j'i 



dij I A Yl a + w, Yl a.j^ Yl 



7=1 



f + J, 



n nK{aj) ^ 

^8,Ba — — H 7^75 r 



EE 



-1 



+ A I + A JJ^ J JJ^ ttj'i 



n 



/ 7^ J, 



A 



Sa,A = 



i=l 



A n +^^11 n ' 



3=1 



A 



i=l 



3=1 / 3=1 

P 



3=1 3=1 / j'T^j 



aj'i 



i=l 



Wi \ A JJ + JJ aji Yl 

3=1 3=1 J 3=1 . 

All the expected values above are obtained numerically. Also, K{aj) — [aj — ^K*{a.j)l\p2\l1, 
with K*{a.j) = [1 — erf (-\/2/Q;j)] exp(2/Q;|), for j — 1, . . . ,p, where erf (•) is the error function 
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given by erf (a;) = (2/^) e Details on erff-) ca n be found in 

(|2007|). For small values of a (|Abramowitz and Stegunl Il970l p. 298) 



Gradshtevn and Rvzhik 



a,- / a? 3q;^ 



For numerical evaluation we recommend the use of (fT2] ) when a < 0.5. 

For A = 0, which corresponds to the independent case, we obtain the Fisher information matrix 

See = nblock-diag{Sc,a, S^^, 2/tt}, 

where Ti^a = 2 diag{af ^, . . . , S^^ = diag{6i, . . . , bp}, with 6^ = [ajK{aj) + 1]/ for 
j = 1, . . . , p. It can be shown that 

\^ee\ 11 ^4^2 ^ 0- 

Therefore, the Fisher information matrix is not singular at A = 0. 

Finally, it is well known that under some mild regularity conditions, the asymptotic behavior 
remains valid if Se is approximated by —Lgg, where —Lgg is the {2p + 1) x (2p + 1) observed 
information matrix evaluated at 6, obtained from Lee. So, in order to avoid numerical integrations, 
one can use —L^g instead of Se to make inference. 
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